home *** CD-ROM | disk | FTP | other *** search
/ Linux Cubed Series 3: Developer Tools / Linux Cubed Series 3 - Developer Tools.iso / devel / lang / forth / pfe-0.000 / pfe-0 / pfe-0.9.13 / src / floating.c < prev    next >
Encoding:
C/C++ Source or Header  |  1995-07-17  |  12.8 KB  |  729 lines
  1. /*
  2.  * This file is part of the portable Forth environment written in ANSI C.
  3.  * Copyright (C) 1995  Dirk Uwe Zoller
  4.  *
  5.  * This library is free software; you can redistribute it and/or
  6.  * modify it under the terms of the GNU Library General Public
  7.  * License as published by the Free Software Foundation; either
  8.  * version 2 of the License, or (at your option) any later version.
  9.  *
  10.  * This library is distributed in the hope that it will be useful,
  11.  * but WITHOUT ANY WARRANTY; without even the implied warranty of
  12.  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
  13.  * See the GNU Library General Public License for more details.
  14.  *
  15.  * You should have received a copy of the GNU Library General Public
  16.  * License along with this library; if not, write to the Free
  17.  * Software Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
  18.  *
  19.  * This file is version 0.9.13 of 17-July-95
  20.  * Check for the latest version of this package via anonymous ftp at
  21.  *    roxi.rz.fht-mannheim.de:/pub/languages/forth/pfe-VERSION.tar.gz
  22.  * or    sunsite.unc.edu:/pub/languages/forth/pfe-VERSION.tar.gz
  23.  * or    ftp.cygnus.com:/pub/forth/pfe-VERSION.tar.gz
  24.  *
  25.  * Please direct any comments via internet to
  26.  *    duz@roxi.rz.fht-mannheim.de.
  27.  * Thank You.
  28.  */
  29. /*
  30.  * floating.c ---     The Optional Floating-Point Word Set
  31.  * (duz 25Jul93)
  32.  */
  33.  
  34. #include "forth.h"
  35. #include "support.h"
  36. #include "compiler.h"
  37. #include "dblsub.h"
  38.  
  39. #include <stdlib.h>
  40. #include <string.h>
  41. #include <ctype.h>
  42. #include <float.h>
  43. #include <math.h>
  44.  
  45. #include "missing.h"
  46.  
  47. Code (d_f_align);
  48.  
  49. #if defined USE_SSCANF        /* define this if you fully trust your scanf */
  50.  
  51. Code (to_float)            /* >FLOAT ( c-addr u --- flag; F: --- r | ) */
  52. /*
  53.  * This is a working solution on most machines.
  54.  * Unfortunately it relies on pretty obscure features of sscanf()
  55.  * which are not truly implemented everywhere.
  56.  */
  57. {
  58.   char *p, buf[80];
  59.   static char *fmt[] =
  60.   {
  61.     "%lf%n %n%d%n$",
  62.     "%lf%*1[DdEe]%n %n%d%n$",
  63.   };
  64.   int i, n, exp, n1, n2, n3;
  65.   double r;
  66.  
  67.   p = (char *) sp[1];
  68.   n = dash_trailing (p, *sp++);
  69.   if (n == 0)
  70.     {
  71.       *--fp = 0;
  72.       *sp = TRUE;
  73.       return;
  74.     }
  75.   store_c_string (p, n, buf, sizeof buf);
  76.   strcat (buf, "$");
  77. #if defined EMX
  78.   /* emx' sscanf(), %lf conversion, doesn't read past 0E accepting the
  79.    * "0" as good number when no exponent follows.  Therefore we change
  80.    * the 'E' to 'D', ugly hack but helps. */
  81.   upper (buf, n);
  82.   if (strchr (buf, 'E'))
  83.     *strchr (buf, 'E') = 'D';
  84. #endif
  85.   if (1 == sscanf (buf, "%lf%n$", &r, &n1) &&
  86.       n == n1)
  87.     {
  88.       *--fp = r;
  89.       *sp = TRUE;
  90.       return;
  91.     }
  92.   for (i = 0; i < DIM (fmt); i++)
  93.     switch (sscanf (buf, fmt[i], &r, &n1, &n2, &exp, &n3))
  94.       {
  95.       case 1:
  96.     if (n < n2)
  97.       break;
  98.     *--fp = r;
  99.     *sp = TRUE;
  100.     return;
  101.       case 2:
  102.     if (n1 != n2 || n < n3)
  103.       break;
  104.     *--fp = r * pow10 (exp);
  105.     *sp = TRUE;
  106.     return;
  107.       }
  108.   *sp = FALSE;
  109. }
  110.  
  111. #else
  112.  
  113. Code (to_float)            /* >FLOAT ( c-addr u --- flag; F: --- r | ) */
  114. /*
  115.  * This is an implementation based on a simple state machine.
  116.  * Uses nothing but simple character manipulation and floating point math.
  117.  */
  118. {
  119.   enum state            /* states of the state machine */
  120.     {
  121.       bpn,            /* before point, maybe sign */
  122.       bp,            /* before point, no more sign (had one) */
  123.       ap,            /* after point */
  124.       exn,            /* exponent, maybe sign */
  125.       ex,            /* exponent, no more sign */
  126.       ts            /* trailing space */
  127.     };
  128.   enum state state = bpn;
  129.   int sign = 1;            /* sign of mantissa */
  130.   long double mant = 0;        /* the mantissa */
  131.   int esign = 1;        /* sign of exponent */
  132.   int exp = 0;            /* the exponent */
  133.   int scale = 0;        /* number of digits after point */
  134.   int n = *sp++;        /* string length */
  135.   char *p = (char *) *sp;    /* points to string */
  136.  
  137.   while (--n >= 0)
  138.     {
  139.       char c = *p++;
  140.  
  141.       switch (state)
  142.     {
  143.     case bpn:
  144.       switch (c)
  145.         {
  146.         case '-':
  147.           sign = -1;
  148.         case '+':
  149.           state = bp;
  150.           continue;
  151.         case '.':
  152.           state = ap;
  153.           continue;
  154.         default:
  155.           if (isspace (c))
  156.         continue;
  157.           if (isdigit (c))
  158.         {
  159.           mant = c - '0';
  160.           state = bp;
  161.           continue;
  162.         }
  163.         }
  164.       goto bad;
  165.     case bp:
  166.       switch (c)
  167.         {
  168.         case '.':
  169.           state = ap;
  170.           continue;
  171.         case '-':
  172.           esign = -1;
  173.         case '+':
  174.           state = ex;
  175.           continue;
  176.         case 'D':
  177.         case 'd':
  178.         case 'E':
  179.         case 'e':
  180.           state = exn;
  181.           continue;
  182.         default:
  183.           if (isspace (c))
  184.         {
  185.           state = ts;
  186.           continue;
  187.         }
  188.           if (isdigit (c))
  189.         {
  190.           mant *= 10;
  191.           mant += c - '0';
  192.           continue;
  193.         }
  194.         }
  195.       goto bad;
  196.     case ap:
  197.       switch (c)
  198.         {
  199.         case '-':
  200.           esign = -1;
  201.         case '+':
  202.           state = ex;
  203.           continue;
  204.         case 'D':
  205.         case 'd':
  206.         case 'E':
  207.         case 'e':
  208.           state = exn;
  209.           continue;
  210.         default:
  211.           if (isspace (c))
  212.         {
  213.           state = ts;
  214.           continue;
  215.         }
  216.           if (isdigit (c))
  217.         {
  218.           mant *= 10;
  219.           mant += c - '0';
  220.           scale--;
  221.           continue;
  222.         }
  223.         }
  224.       goto bad;
  225.     case exn:
  226.       switch (c)
  227.         {
  228.         case '-':
  229.           esign = -1;
  230.         case '+':
  231.           state = ex;
  232.           continue;
  233.         default:
  234.           if (isspace (c))
  235.         {
  236.           state = ts;
  237.           continue;
  238.         }
  239.           if (isdigit (c))
  240.         {
  241.           exp = c - '0';
  242.           state = ex;
  243.           continue;
  244.         }
  245.         }
  246.       goto bad;
  247.     case ex:
  248.       if (isspace (c))
  249.         {
  250.           state = ts;
  251.           continue;
  252.         }
  253.       if (isdigit (c))
  254.         {
  255.           exp *= 10;
  256.           exp += c - '0';
  257.           continue;
  258.         }
  259.       goto bad;
  260.     case ts:
  261.       if (isspace (c))
  262.         continue;
  263.       goto bad;
  264.     }
  265.     }
  266.   *--fp = sign * mant * pow10 (scale + esign * exp);
  267.   *sp = TRUE;
  268.   return;
  269. bad:
  270.   *sp = FALSE;
  271.   return;
  272. }
  273.  
  274. #endif
  275.  
  276. Code (d_to_f)
  277. {
  278.   int sign;
  279.   double res;
  280.  
  281.   if (sp[0] < 0)
  282.     sign = 1, dnegate ((dCell *) &sp[0]);
  283.   else
  284.     sign = 0;
  285. #if Linux
  286.   /* slackware 2.2.0.1 (at least) has a bug in ldexp()  */
  287.   res = (uCell) sp[0] * ((double)(1<<31) * 2) + (uCell) sp[1];
  288. #else
  289.   res = ldexp ((uCell) sp[0], CELLBITS) + (uCell) sp[1];
  290. #endif
  291.   sp += 2;
  292.   *--fp = sign ? -res : res;
  293. }
  294.  
  295. Code (f_store)
  296. {
  297.   *(double *) *sp++ = *fp++;
  298. }
  299.  
  300. Code (f_star)
  301. {
  302.   fp[1] *= fp[0];
  303.   fp++;
  304. }
  305.  
  306. Code (f_plus)
  307. {
  308.   fp[1] += fp[0];
  309.   fp++;
  310. }
  311.  
  312. Code (f_minus)
  313. {
  314.   fp[1] -= fp[0];
  315.   fp++;
  316. }
  317.  
  318. Code (f_slash)
  319. {
  320.   fp[1] /= fp[0];
  321.   fp++;
  322. }
  323.  
  324. Code (f_zero_less)
  325. {
  326.   *--sp = FLAG (*fp++ < 0);
  327. }
  328.  
  329. Code (f_zero_equal)
  330. {
  331.   *--sp = FLAG (*fp++ == 0);
  332. }
  333.  
  334. Code (f_less_than)
  335. {
  336.   *--sp = FLAG (fp[1] < fp[0]);
  337.   fp += 2;
  338. }
  339.  
  340. Code (f_to_d)
  341. {
  342.   double a, hi, lo;
  343.   int sign;
  344.  
  345.   if ((a = *fp++) < 0)
  346.     sign = 1, a = -a;
  347.   else
  348.     sign = 0;
  349.   lo = modf (ldexp (a, -CELLBITS), &hi);
  350.   sp -= 2;
  351.   sp[0] = (uCell) hi;
  352.   sp[1] = (uCell) ldexp (lo, CELLBITS);
  353.   if (sign)
  354.     dnegate ((dCell *) &sp[0]);
  355. }
  356.  
  357. Code (f_fetch)
  358. {
  359.   *--fp = *(double *) *sp++;
  360. }
  361.  
  362. void
  363. f_constant_runtime (void)
  364. {
  365.   *--fp = *(double *) dfaligned ((Cell) PFA);
  366. }
  367.  
  368. Code (f_constant)
  369. {
  370.   header (f_constant_runtime, 0);
  371.   d_f_align_ ();
  372.   FCOMMA (*fp++);
  373. }
  374.  
  375. Code (f_depth)
  376. {
  377.   *--sp = sys.f0 - fp;
  378. }
  379.  
  380. Code (f_drop)
  381. {
  382.   fp++;
  383. }
  384.  
  385. Code (f_dup)
  386. {
  387.   fp--;
  388.   fp[0] = fp[1];
  389. }
  390.  
  391. code (f_literal_execution)
  392. {
  393.   POP (double, ip, *--fp);
  394. }
  395.  
  396. code (f_literal)
  397. {
  398.   if (STATE)
  399.     {
  400. #if DFLOAT_ALIGN > CELL_ALIGN
  401.       if (DFALIGNED (DP))
  402.     compile2 ();
  403. #endif
  404.       compile1 ();
  405.       FCOMMA (*fp++);
  406.     }
  407. }
  408. COMPILES2 (f_literal, f_literal_execution, noop,
  409.        SKIPS_FLOAT, DEFAULT_STYLE);
  410.  
  411. Code (floor)
  412. {
  413.   *fp = floor (*fp);
  414. }
  415.  
  416. Code (f_max)
  417. {
  418.   if (fp[0] > fp[1])
  419.     fp[1] = fp[0];
  420.   fp++;
  421. }
  422.  
  423. Code (f_min)
  424. {
  425.   if (fp[0] < fp[1])
  426.     fp[1] = fp[0];
  427.   fp++;
  428. }
  429.  
  430. Code (f_negate)
  431. {
  432.   *fp = -*fp;
  433. }
  434.  
  435. Code (f_over)
  436. {
  437.   fp--;
  438.   fp[0] = fp[2];
  439. }
  440.  
  441. Code (f_rot)
  442. {
  443.   double h = fp[2];
  444.  
  445.   fp[2] = fp[1];
  446.   fp[1] = fp[0];
  447.   fp[0] = h;
  448. }
  449.  
  450. Code (f_round)
  451. {
  452.   *fp = floor (*fp + 0.5);
  453. }
  454.  
  455. Code (f_swap)
  456. {
  457.   double h = fp[1];
  458.  
  459.   fp[1] = fp[0];
  460.   fp[0] = h;
  461. }
  462.  
  463. void
  464. f_variable_runtime (void)
  465. {
  466.   *--sp = dfaligned ((Cell) PFA);
  467. }
  468.  
  469. Code (f_variable)
  470. {
  471.   header (f_variable_runtime, 0);
  472.   d_f_align_ ();
  473.   FCOMMA (0.);
  474. }
  475.  
  476. Code (represent)        /* with help from Lennart Benshop */
  477. {
  478.   char *p, buf[0x80];
  479.   int u, log, sign;
  480.   double f;
  481.  
  482.   f = *fp++;
  483.   p = (char *) sp[1];
  484.   u = sp[0];
  485.   sp--;
  486.  
  487.   if (f < 0)
  488.     sign = TRUE, f = -f;
  489.   else
  490.     sign = FALSE;
  491.   if (f != 0)
  492.     {
  493.       log = (int) floor (log10 (f)) + 1;
  494.       f *= pow10 (-log);
  495.       if (f + 0.5 * pow10 (-u) >= 1)
  496.     f /= 10, log++;
  497.     }
  498.   else
  499.     log = 0;
  500.   sprintf (buf, "%0.*f", u, f);
  501.   memcpy (p, buf + 2, u);
  502.  
  503.   sp[2] = log;
  504.   sp[1] = sign;
  505.   sp[0] = TRUE;
  506. }
  507.  
  508. /************************************************************************/
  509. /* Floating point extension words:                                      */
  510. /************************************************************************/
  511.  
  512. Code (d_f_align)
  513. {
  514.   while (!DFALIGNED (DP))
  515.     *DP++ = 0;
  516. }
  517.  
  518. Code (d_f_aligned)
  519. {
  520.   sp[0] = dfaligned (sp[0]);
  521. }
  522.  
  523. Code (d_float_plus)
  524. {
  525.   *sp += sizeof (double);
  526. }
  527.  
  528. Code (d_floats)
  529. {
  530.   *sp *= sizeof (double);
  531. }
  532.  
  533. Code (f_star_star)
  534. {
  535.   fp[1] = pow (fp[1], fp[0]);
  536.   fp++;
  537. }
  538.  
  539. Code (f_dot)
  540. {
  541.   outf ("%.*f ", PRECISION, *fp++);
  542. }
  543.  
  544. Code (f_abs)
  545. {
  546.   if (*fp < 0)
  547.     *fp = -*fp;
  548. }
  549.  
  550. Code (f_e_dot)            /* with help from Lennart Benshop */
  551. {
  552.   double f = fabs (*fp);
  553.   double h = 0.5 * pow10 (-PRECISION);
  554.   int n;
  555.  
  556.   if (f == 0)
  557.     n = 0;
  558.   else if (f < 1)
  559.     {
  560.       h = 1 - h;
  561.       for (n = 3; f * pow10 (n) < h; n += 3);
  562.     }
  563.   else
  564.     {
  565.       h = 1000 - h;
  566.       for (n = 0; h <= f * pow10 (n); n -= 3);
  567.     }
  568.   outf ("%+*.*fE%+03d ", PRECISION + 5, PRECISION,
  569.     *fp++ * pow10 (n), -n);
  570. }
  571.  
  572. Code (f_s_dot)
  573. {
  574.   outf ("%.*E ", PRECISION, *fp++);
  575. }
  576.  
  577. Code (f_proximate)
  578. {
  579.   double a, b, c;
  580.  
  581.   a = fp[2];
  582.   b = fp[1];
  583.   c = fp[0];
  584.   fp += 3;
  585. #if 0
  586.   sp--;
  587.   if (c > 0)
  588.     *sp = FLAG (fabs (a - b) < c);
  589.   else if (c < 0)
  590.     *sp = FLAG (fabs (a - b) < -c * (fabs (a) + fabs (b)));
  591.   else
  592.     *sp = FLAG (memcmp (&a, &b, sizeof (double)) == 0);
  593.  
  594. #else
  595.   *--sp = FLAG
  596.     (c > 0 ? fabs (a - b) < c :
  597.      c < 0 ? fabs (a - b) < -c * (fabs (a) + fabs (b))
  598.      : a == b);
  599. #endif
  600. }
  601.  
  602. Code (set_precision)
  603. {
  604.   PRECISION = *sp++;
  605. }
  606.  
  607. Code (s_f_store)
  608. {
  609.   *(float *) *sp++ = *fp++;
  610. }
  611.  
  612. Code (s_f_fetch)
  613. {
  614.   *--fp = *(float *) *sp++;
  615. }
  616.  
  617. Code (s_float_plus)
  618. {
  619.   *sp += sizeof (float);
  620. }
  621.  
  622. Code (s_floats)
  623. {
  624.   *sp *= sizeof (float);
  625. }
  626.  
  627. /* simple mappings to the ANSI-C library: */
  628. /* *INDENT-OFF* */
  629. Code (f_acos)    { *fp = acos (*fp); }
  630. Code (f_acosh)    { *fp = acosh (*fp); }
  631. Code (f_alog)    { *fp = pow10 (*fp); }
  632. Code (f_asin)    { *fp = asin (*fp); }
  633. Code (f_asinh)    { *fp = asinh (*fp); }
  634. Code (f_atan)    { *fp = atan (*fp); }
  635. Code (f_atan2)    { fp [1] = atan2 (fp [1], fp [0]); fp++; }
  636. Code (f_atanh)    { *fp = atanh (*fp); }
  637. Code (f_cos)    { *fp = cos (*fp); }
  638. Code (f_cosh)    { *fp = cosh (*fp); }
  639. Code (f_exp)    { *fp = exp (*fp); }
  640. Code (f_expm1)    { *fp = exp (*fp) - 1.0; }
  641. Code (f_ln)    { *fp = log (*fp); }
  642. Code (f_lnp1)    { *fp = log (*fp + 1.0); }
  643. Code (f_log)    { *fp = log10 (*fp); }
  644. Code (f_sin)    { *fp = sin (*fp); }
  645. Code (f_sincos)    { --fp; fp [0] = cos (fp [1]); fp [1] = sin (fp [1]); }
  646. Code (f_sinh)    { *fp = sinh (*fp); }
  647. Code (f_sqrt)    { *fp = sqrt (*fp); }
  648. Code (f_tan)    { *fp = tan (*fp); }
  649. Code (f_tanh)    { *fp = tanh (*fp); }
  650.  
  651.  
  652. LISTWORDS (floating) =
  653. {
  654.   CO (">FLOAT",        to_float),
  655.   CO ("D>F",        d_to_f),
  656.   CO ("F!",        f_store),
  657.   CO ("F*",        f_star),
  658.   CO ("F+",        f_plus),
  659.   CO ("F-",        f_minus),
  660.   CO ("F/",        f_slash),
  661.   CO ("F0<",        f_zero_less),
  662.   CO ("F0=",        f_zero_equal),
  663.   CO ("F<",        f_less_than),
  664.   CO ("F>D",        f_to_d),
  665.   CO ("F@",        f_fetch),
  666.   CO ("FALIGN",        d_f_align),
  667.   CO ("FALIGNED",    d_f_aligned),
  668.   CO ("FCONSTANT",    f_constant),
  669.   CO ("FDEPTH",        f_depth),
  670.   CO ("FDROP",        f_drop),
  671.   CO ("FDUP",        f_dup),
  672.   CS ("FLITERAL",    f_literal),
  673.   CO ("FLOAT+",        d_float_plus),
  674.   CO ("FLOATS",        d_floats),
  675.   CO ("FLOOR",        floor),
  676.   CO ("FMAX",        f_max),
  677.   CO ("FMIN",        f_min),
  678.   CO ("FNEGATE",    f_negate),
  679.   CO ("FOVER",        f_over),
  680.   CO ("FROT",        f_rot),
  681.   CO ("FROUND",        f_round),
  682.   CO ("FSWAP",        f_swap),
  683.   CO ("FVARIABLE",    f_variable),
  684.   CO ("REPRESENT",    represent),
  685.   /* floating point extension words */
  686.   CO ("DF!",        f_store),
  687.   CO ("DF@",        f_fetch),
  688.   CO ("DFALIGN",    d_f_align),
  689.   CO ("DFALIGNED",    d_f_aligned),
  690.   CO ("DFLOAT+",    d_float_plus),
  691.   CO ("DFLOATS",    d_floats),
  692.   CO ("F**",        f_star_star),
  693.   CO ("F.",        f_dot),
  694.   CO ("FABS",        f_abs),
  695.   CO ("FACOS",        f_acos),
  696.   CO ("FACOSH",        f_acosh),
  697.   CO ("FALOG",        f_alog),
  698.   CO ("FASIN",        f_asin),
  699.   CO ("FASINH",        f_asinh),
  700.   CO ("FATAN",        f_atan),
  701.   CO ("FATAN2",        f_atan2),
  702.   CO ("FATANH",        f_atanh),
  703.   CO ("FCOS",        f_cos),
  704.   CO ("FCOSH",        f_cosh),
  705.   CO ("FE.",        f_e_dot),
  706.   CO ("FEXP",        f_exp),
  707.   CO ("FEXPM1",        f_expm1),
  708.   CO ("FLN",        f_ln),
  709.   CO ("FLNP1",        f_lnp1),
  710.   CO ("FLOG",        f_log),
  711.   CO ("FS.",        f_s_dot),
  712.   CO ("FSIN",        f_sin),
  713.   CO ("FSINCOS",    f_sincos),
  714.   CO ("FSINH",        f_sinh),
  715.   CO ("FSQRT",        f_sqrt),
  716.   CO ("FTAN",        f_tan),
  717.   CO ("FTANH",        f_tanh),
  718.   CO ("F~",        f_proximate),
  719.   SC ("PRECISION",    PRECISION),
  720.   CO ("SET-PRECISION",    set_precision),
  721.   CO ("SF!",        s_f_store),
  722.   CO ("SF@",        s_f_fetch),
  723.   CO ("SFALIGN",    align),
  724.   CO ("SFALIGNED",    aligned),
  725.   CO ("SFLOAT+",    s_float_plus),
  726.   CO ("SFLOATS",    s_floats),
  727. };
  728. COUNTWORDS (floating, "Floating point + extensions");
  729.